%-------------------------------------------------------------------------------
% $Id: wav1d_finite_diff.m,v 1.4 2011/06/25 17:34:43 paul Exp $
% $Date: 2011/06/25 17:34:43 $
% $Author: paul $
%
%% Explicit Finite Difference Solution to 1D wave equation of a vibrating string
%
% du{tt} = (c^2)*du{xx}
% Initial Conditions: u(x,0) = f(x)
%                     du(x,0){t} = g(x)
% Boundary Conditions: u(0,t) = u(L,t) = 0
%-------------------------------------------------------------------------------
function wav1d_finite_diff(arg_M, arg_N)
% Explicit Finite Difference Solution to 1D wave equation

if nargin < 2
  M = 128;                       % number of samples on x-axis
  N = 256;                       % number of samples in time interval T
  msg = sprintf('INFO:\tUsing the default arguments M=%d N=%d', M,N);
  disp(msg);
else           
  M = arg_M;
  N = arg_N;
  msg = sprintf('INFO:\tUser arguments M=%d N=%d"', M,N);
  disp(msg);
end

T = 2;                         % time interval (0,T)
c = 1.0;                       % constant tension
length = 1.0;                  % length of string on x-axis
h = length / M;                % dx
k = T / N;                     % dt
lamda = c * k/h;               % lamda =< 1 for stability

msg = sprintf('INFO:\tParameters for Finite Difference Solution:');
disp(msg); 
msg = sprintf('INFO:\tT=%f c=%f h=%f k=%f lamda=%f',T,c,h,k,lamda);
disp(msg); 
  
prompt = sprintf('INPUT:\tPress return to continue');
response=input(prompt);

if (lamda > 1)
  warning = sprintf('WARNING:lamda = %f. PDE solution unstable if lamda > 1', lamda);
  disp(warning);
  hint = sprintf('INFO:\tAdjust parameters M, N, c and T for convergence');
  disp(hint);
  prompt = sprintf('INPUT:\tPress return to continue');
  response=input(prompt);
end

%-------------------------------------------------------------------------------
% Setup AVI file and parameters
aviobj = avifile('wav1d_finite_diffs.avi');
aviobj.compression = 'None';   % no other option for Unix
aviobj.fps = 5;                % frames per second
aviobj.quality = 10;           % low quality is OK

%-------------------------------------------------------------------------------
% Animation loop
fig = figure;

lamda_sq = lamda ^ 2;          % precompute some values outside of loop!
factor = 1 - lamda_sq;           
factor2 = 2*factor;

for j=1:N
  if (j==1)
    % initial conditions of y
    x = 0:h:1;
    y = zeros(1,M+1);              % extra padding element at end
    yd1 = zeros(1,M+1);            % last value of y 
    yd2 = zeros(1,M+1);            % one earlier value of y
    
    M1 = floor(0.2 * M);           % window of (0,0.2) for initial displacement
                                   % round M1 to nearest integer
    
    y(2:M1) = exp(-1./(1-(10*x(2:M1)-1).^2));
    g = zeros(1,M+1);              % zero initial velocity 
    
    % boundary conditions - maintain fixed end points y(1) = y(M) = 0
  elseif (j==2)
    % first iteration
    y(2:M-1) = factor * yd1(2:M-1) + (lamda_sq/2) * (yd1(3:M)+yd1(1:M-2)) + k*g(2:M-1);
    y(1) = 0; y(M) = 0;    
  else
    % remaining iterations
    y(2:M-1) = factor2 * yd1(2:M-1) + lamda_sq * (yd1(3:M)+yd1(1:M-2)) - yd2(2:M-1);    
    y(1) = 0; y(M) = 0;
  end
  
  % update delayed versions of string displacement state (amplitude) y 
  yd2 = yd1;
  yd1 = y;

  plot([1:M],y(1:M));
  axis([1 M -0.4 0.4]);
  title('1D Wave Equation Finite Differences Solution');
  frame = getframe(fig); 
  aviobj = addframe(aviobj,frame);
end 

%-------------------------------------------------------------------------------
% Clean up rendering frame and AVI file
close(fig);
aviobj = close(aviobj);